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ABSTRACT. We present a method for speeding up numerical calculations of a light curve for a stellar 
occultation by a planetary atmosphere with an arbitrary atmospheric model that has spherical symmetry. 

This improved speed makes least-squares fitting for model parameters practical. Our method takes as input 
several sets of values for the first two radial derivatives of the refractivity at different values of model 
parameters, and interpolates to obtain the light curve at intermediate values of one or more model 
parameters. It was developed for small occulting bodies such as Pluto and Triton, but is applicable to planets 
of all sizes. We also present the results of a series of tests showing that our method calculates light curves 
that are correct to an accuracy of 10 4 of the unocculted stellar flux. The test benchmarks are (i) an 
atmosphere with a Mr dependence of temperature, which yields an analytic solution for the light curve, (ii) 
an atmosphere that produces an exponential refraction angle, and (iii) a small-planet isothermal model. With 
our method, least-squares fits to noiseless data also converge to values of parameters with fractional errors 
of no more than 10 4 , with the largest errors occurring in small planets. These errors are well below the 
precision of the best stellar occultation data available. Fits to noisy data had formal errors consistent with 
the level of synthetic noise added to the light curve. We conclude: (i) one should interpolate refractivity 
derivatives and then form light curves from the interpolated values, rather than interpolating the light curves 
themselves; (ii) for the most accuracy, one must specify the atmospheric model for radii many scale heights 
above half light; and (iii) for atmospheres with smoothly varying refractivity with altitude, light curves can 
be sampled as coarsely as two points per scale height. 



1. INTRODUCTION 

Since Baum and Code’s observations of a stellar occul- 
tation by Jupiter in 1952 (Baum and Code 1953), the analy- 
sis of occultation data to learn about the structure of plan- 
etary atmospheres has developed along two lines: numerical 
inversion (Kovalevsky and Link 1969; Wasserman and 
Veverka 1973; French et al. 1978) and model fitting (Baum 
and Code 1953; Elliot and Young 1992). In principle, for 
atmospheres in hydrostatic equilibrium, the inversion ap- 
proach will yield the temperature profile, as well as pressure 
and number density profiles (if one knows the planetary 
gravity and mean molecular weight — which effectively 
means knowing the composition). However, small uncer- 
tainties in the initial conditions cause significant errors in 
the results until the inversion calculation has been carried 
out for several scale heights (Wasserman and Veverka 
1973; French et al. 1978; Roques et al. 1994). By then, 
however, the occultation light-curve flux is low and errors 
in the zero flux level can significantly affect the results. 
Noise in the light curve also adds spurious structure to the 
inversion results. In spite of these difficulties, inversion is 
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the best method of analysis when one does not have an a 
priori model of the atmospheric structure. 

The modeling approach has been limited to isothermal 
atmospheres (Baum and Code 1953) and atmospheres with 
power-law thermal gradients (Elliot and Young 1992, here- 
after EY92). However, our theoretical understanding of at- 
mospheric structure has progressed far beyond just estimat- 
ing a scale height, and the thermal gradients in atmospheres 
would not necessarily follow a power law. For example, 
physical models of Triton’s atmosphere (Krasnopolsky et al. 
1993; Strobel and Summers 1995) include the physical ef- 
fects of electron heating, solar radiative heating of CH 4 , 
thermal conduction, and radiative cooling by CH 4 and CO. 
With a stellar occultation data set of high quality (e.g., 
Olkin 1996), one can hope to test the validity of these 
models, and — for those models found to be valid — to use 
the data to estimate values of atmospheric parameters, such 
as the CH 4 and CO fractions in the atmosphere. 

To carry out these analyses, one would like to fit an 
appropriate atmospheric model to the occultation light curve 
by least squares, in order to find the most likely values of 
the model parameters and their errors. Because of the com- 
plexity of the atmospheric models, a numerical approach is 
required: for each set of model parameters, a model light 
curve must be calculated, and for each parameter being fit, 
one must calculate a numerical derivative for each point on 
the light curve by incrementing the parameter and calculat- 
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ing another model light curve. Each time the physical 
model is recalculated (e.g., Strobel and Summers 1995), this 
process involves equations that can only be solved by nu- 
merical iteration to find self-consistent values, which, at 
present, takes several hours (Strobel, personal communica- 
tion). Hence computing a light curve and its numerical 
derivatives from such a model for each parameter fit can 
involve an enormous amount of computing time for each set 
of parameters. This computing time can make least-squares 
fitting for the best parameter values cumbersome at best, 
and often impractical. 

In this paper we present a method for speeding up the 
process of calculating light curves from complicated atmo- 
spheric models by interpolating a grid of precalculated 
models. In calculating the light curve, we include geometri- 
cal effects appropriate to small planets (EY92). Our model 
is thus applicable to occulting bodies like Pluto and Triton 
as well as to larger planets. We also allow for the effects of 
extinction, though we do not present results of tests that 
include it. We start with one or more profiles of values for 
the first two radial derivatives of the refractivity at radii 
spanning the region probed by the occultation, and we in- 
terpolate between these values to calculate a model light 
curve for some set of parameters. We present a method for 
interpolating between points along a single refractivity pro- 
file to make a model light curve from a single set of pa- 
rameters and for interpolating between models with differ- 
ing values of parameters to make light curves at arbitrary 
intermediate values of the model parameters. 

In deriving the equations for converting a refractivity 
profile into a light curve, we closely follow Sec. 2 of EY92. 
As they did, we use the geometrical optics approximation 
and make some simplifying assumptions: (i) the light re- 
ceived at any time comes from only one point on the plan- 
etary limb, (ii) the atmosphere is spherically symmetric. 
These assumptions are not critical to our method, however, 
and we make them only to simplify the numerical compu- 
tations presented here. 

2. THE MODEL 


2.1 Analytic Basis 

We assume that parallel, monochromatic light lays are 
incident on a spherically symmetric planetary atmosphere 
from the left in Fig. 1. The observer’s plane is perpendicu- 
lar to the path of the light rays. In this plane, y is a radial 
coordinate measuring the observer’s position relative to the 
point through which a line connecting the occulted star and 
the center of the planet would pass. EY92 used a coordinate 
p = | y j. Using >• instead of p enables us to more easily 
include contributions to the flux from both the near and far 
limbs (although in this paper we shall be working with only 
the near-limb flux). We refer to y as a function of r, which 
marks the position of arrival in the observer’s plane of a 
light ray with a radius of the closest approach to the planet 
of r. We also refer to y as a function of t, where y(f) 
marks the position of the observer in the plane as a function 
of time. This function can be determined from the geocen- 
tric planetary ephemeris and the motion of the observer 



Fig. 1 — Stellar occultation by a planetary atmosphere. Starlight encounters 
a planetary atmosphere and is bent by the refractivity gradient in the 
atmosphere. Since the refraction increases exponentially with depth into 
the atmosphere, refraction causes two neighboring rays to be bent by a 
slightly different angle, and this differential separation increases by an 
amount proportional to distance from the planet. This effect causes the star 
to dim as seen by a distant observer (after EY92). 


relative to the center of the Earth. The refraction angle, 
6{r), is defined in the usual sense to be positive above the 
axis, so refraction by the atmosphere produces negative 
values. 

With these assumptions, the flux from the occulted star 
can be broken down into three terms: (i) the differential 
refraction of the light rays at different depths in the atmo- 
sphere, (ii) the partial focusing of the light rays in the plane 
perpendicular to the path of the ray, and (iii) the extinction 
of the light due to absorption in the atmosphere. From 
EY92, this yields the following equation for the flux, £, as 
a function of r : 


Sir) = 


1 l 

| J + D[dO(r)fdr]\ |l + D0(r)/r 

differential refraction focusing 


exp[- r obs (r)] 
extinction 

(i) 


where r obs (r) is the optical depth along the path of the ray, 
D is the distance between the occulting body and the ob- 
server, and 0(r) is the angle through which a light ray is 
refracted as a function of the closest approach, r, to the 
center of the occulting body. If we made no assumptions 
about the refractivity in the atmosphere, we would need to 
include a factor of one over the index of refraction in the 
calculation of 6{r) (Born and Wolf 1964). To simplify the 
calculation we assume that the refractivity 
v <§ l and that 1 6>| <§ 1 throughout the atmosphere. Under 
this assumption, the following approximations are valid. 


6(r) 


' x r dv(r') 

— r~ . , dx 
— oc r dr 


(2) 


dO(r) 

dr 


' * [ x 2 dv{r') 
_x (r') 3 dr' 


r 2 d 2 v(r') | 


(3) 


where x is a coordinate along the path of the ray measured 
from its point of closest approach and r f is the radial 
distance from the planet along the path of the ray (see Fig. 
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1). The variables r, r\ and x satisfy the relation 
r fl = r 2 + x 2 . 

In this paper we shall be referring to (p cy] (>’), the unfo- 
cused stellar flux from one limb that is not affected by 
extinction. The subscript “cyl” reminds us that this is the 
flux that would be seen refracted by a cylindrical atmo- 
sphere (one with no focusing in the plane perpendicular to 
the ray path). The range of <f> cy \(y) corresponds to a nor- 
malized stellar flux, which equals 1.0 when the starlight is 
unaffected by the atmosphere and 0.0 when the stellar flux 
is totally occulted. Since focusing and extinction are not 
included, our equation for <£ cyl (y) includes only the first 
term of Eq. (1): 


4*cy\(y ) 


1 +D 


dO(r)\ 

dr 


where y is related to r by 

y = r + Dd{r). 


(4) 


(5) 


The derivatives of the refractivity required to calculate 6{r) 
and d0(r)ldr are supplied by the atmospheric model. Note 
that only the first and second derivatives of the refractivity 
are needed to calculate the light curve, and not the refrac- 
tivity itself. 

For atmospheres with v{r ) monotonically decreasing 
with r, the derivative in Eqs. (1) and (4) will always be 
positive, and <f> cy[ will be finite. However, this condition 
does not hold for all atmospheric models (and all real at- 
mospheres), so that <f) cyl can have infinite values when the 
denominator in Eq. (4) equals zero. We have not addressed 
this problem in the method presented here. This problem 
could be addressed in the future by integrating <f> cy \(y) over 
a small interval Ay. This integral will be finite since the 
total light flux is finite. Also, the proposed integration pro- 
cedure for treating points where Eq. (4) becomes infinite 
would correspond to a real light curve, since any observed 
light curve is always integrated (by the data recording 
equipment) over some finite interval of distance within the 
occultation shadow. 


2.2 Numerical Implementation 


The next step is to evaluate the equations developed in 
the previous section using only a set of values of the first 
two derivatives of v(r r ) at discrete radii. In order to use 
Eq. (1) to calculate the flux, we need to put dv{r’)!dr 9 
and d 2 v(r' )/dr' 2 into integrable forms. We do this by 
interpolating between the values we have for these quanti- 
ties. The interpolation functions are only valid out to .v max , 
where 


x 2 - r 2 - 

max ' max 


( 6 ) 


Hence, we approximate 0 and dOfdr by numerically 
integrating Eqs. (2) and (3) out to ±x max . Since the angle 
of bending in atmospheres falls off exponentially with in- 
creasing r, in principle we can choose x max so that the 
contribution to the integral beyond it is as small as desired. 


To calculate a model light curve that we can compare 
with occultation data we need the flux, and hence 9 and 
d 6! dr, as a function of time. Since we have 9 and d9ldr 
as functions of radius, we want to find the radius as a 
function of time. We do this in two steps. First we choose 
a set of values of r (we used equally spaced points for 
simplicity) and use our formula for 9 in Eq. (5) to make a 
list of (r,y) pairs. We use these pairs to make an interpo- 
lation function for r(y ). Next we determine the function 
y(t) from the details of the occultation we observed. For 
small planet occultations, EY92 ? s Eq. (5.1) is often appro- 
priate. For larger planets and/or longer occultations, the 
linear approximation for y(r) may not be adequate, and we 
must calculate it from the planetary ephemeris and the mo- 
tion of the observer for each desired time, t (Elliot et al. 
1993). 

Once we have r(y) and y(t), we are able to calculate 
the radius corresponding to each time in which we are 
interested. We can then use our equations for 9 and dBfdr 
to calculate the first two terms in Eq. (1) at these radii, 
which is the flux we are interested in (see Fig. 2). 

2.3 Interpolations 

An alternative method for calculating light curves for 
parameter values between those of the given refractivity 
derivative profiles would be to calculate the light curve for 
each refractivity profile, then interpolate between the light 
curves themselves, rather than interpolating the refractivi- 
ties. However, interpolation assumes that the function being 
approximated is smooth; when there are discontinuities or 
sharp changes in the function, interpolation can introduce 
large errors. Light curves can have discontinuities when, for 
example, the surface of the planet moves in front of the 
occulted star. Geometrical effects near the center of the 
shadow can also cause abrupt changes in the light curve, 
such as the “central flash” (Elliot et al. 1977). Since the 
derivatives of the refractivity change more smoothly, the 
better course is to interpolate these derivatives and then 
calculate the light curve (rather than interpolating a grid of 
light curves). Another reason for interpolating the deriva- 
tives of the refractivity is that if one has data from several 
stations for the same occultation by a planet with a spheri- 
cally symmetric atmosphere, one need calculate only one 
interpolation function for the refractivity derivatives (since 
the refractivity as a function of planetary radius will be the 
same for all stations). If one interpolated the light curves 
instead, different interpolation functions would be needed 
for each station, since the light curves will look very dif- 
ferent at stations that probed the atmosphere at different 
minimum distances from the center of the occultation 
shadow. 

The interpolation functions we made for the derivatives 
of the refractivity profiles needed to be multidimensional. 
Many methods for making such interpolation functions re- 
quire a rectangular grid of points as inputs. Rather than 
requiring that each profile have values for the derivatives at 
the same set of radii, we chose to make the interpolation 
functions in two steps. First, we made one-dimensional in- 
terpolation functions for the derivatives from each profile. 
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from Eq. (2.1) 


limb points 


Fig. 2 — Calculation flow chart. This diagram summarizes the steps in our method for numerically calculating the light curve from refractivity derivative 
profiles. The arrows represent procedures and the solid boxes represent the input to these procedures and their results. Dotted lines indicate optional steps. 


We then used these interpolation functions to calculate the 
derivatives of the refractivity for each set of parameter 
values at a common set of radii, thereby giving us a rect- 
angular grid of points that we could use to make our mul- 
tidimensional interpolation functions. 

2.4 Boundary Conditions at Large Radii 

It is often necessary to fit a long section of the light 
curve when the star is not occulted in addition to the oc- 
cultation event itself in order to establish the full scale 
brightness of the star. With the method described above, 
fitting a long upper baseline requires that values of the 
refractivity be calculated out to very large radii. If calculat- 
ing such values is time consuming, we can make a further 
approximation for the flux beyond r max to extend the 
method with less computation time. We reject the simplest 
possible approximation of setting the flux to 1.0 beyond 
r max because this introduces a discontinuity in the flux, 
which can interfere with least-squares fitting of these mod- 
els to occultation data. Instead, we approximate the values 
at the ends using the model described in Secs. 3 and 4 of 
EY92, calculating their model to only first order and assum- 
ing that the outer atmosphere is isothermal. A detailed ex- 
planation of how to do this is presented by Chamberlain 
(1996). 

2.5 Software Implementation 

We chose to implement our method as a package in 
Mathematica™ 2.2 (Wolfram 1991) on a Power Macintosh 
8100/80 running a remote kernel on an HP9000 735/125. 
Mathematica™ compromised some computing speed for 
ease of development. The speed penalty could be amelio- 
rated by implementing our method either fully in the C 


programming language, or in something like Mathematica’ s 
MathLink™ protocol, which allows calls to C programs 
from Mathematica™ notebooks. From tests comparing the 
speed of C programs to that of Mathematica™ notebooks, 
we estimate that our method would run about 100 times 
faster on the same machine if we implemented it in C rather 
than in Mathematica™. 

To increase speed as much as possible within Math- 
ematica™, we used Mathematica’ s internal numerical inte- 
gration and interpolation functions. Mathematica™ uses 
piece-wise continuous polynomials to interpolate and adap- 
tive quadrature to numerically integrate. Table I details 
what we interpolated and what order polynomials we used. 


3. TESTS 

We performed several tests of our method. We tested its 
accuracy for calculating single light curves, and then for 
multiple light curves (as would be needed for least-squares 
fitting of atmospheric models). Also we note the computing 
time for a typical calculation. 


3.1 Single Light Curves 

The three model atmospheres used in our accuracy tests 
were: (i) an analytic model atmosphere, for which we can 
specify a power-law temperature profile that corresponds to 
a light curve that can be written with familiar functions 
(Eshleman and Gurrola 1993); (ii) a tightly bound atmo- 
sphere, close to that originally given by Baum and Code 
(1953); (iii) a loosely bound atmosphere from the small- 
planet model of EY92. 
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Table 1 

Summary of Interpolated Functions 


Function 

Order of 

Interpolated 

Spacing of Points 

Comments 

interpolated 

interpolation 

variable 

(per scale height) 


dv(r,{Pi})/dr, 
d 2 v(r y {p i })/dr 2 

3 

r 

2-10 

Interpolated the natural log of the function 
Included derivatives for dv(r)!dr 

dv{r y {p t })idr. 

3 

r. 

30-40 

Interpolated the natural log of the function 

d 2 v(r,{Pi})/dr 2 


{p,Y 

b 

Spacing of points varied 

r{y) 

3 

y 

50 

Spacing of points varied 


a {/?,} represents one or more parameters of the atmospheric model. 

b Optimum spacing of atmospheric parameter values is highly dependent on the parameter. 


3.1.1 Analytic Atmosphere 

To describe the degree to which a planetary atmosphere 
is bound, we use the parameter \ h9 which, as defined by 
Eq. (3.9) of EY92, is equal to the ratio of the gravitational 
potential to kT at the half-light radius, r h :X h = GM/xm/ 
kTr h , where G is the gravitational constant, M is the mass 
of the planet, fi is the molecular weight, m is the atomic 
mass unit, k is Boltzmann’s constant, and T is the tempera- 
ture. We made a refractivity profile using the model speci- 
fied in the Appendix (A.l) with X h = 100 and r h = 1500 
km and compared the light curve we calculated with our 
method to a light curve calculated analytically with the 
model. The light curve and the differences between our 
method and the model light curve are shown in Figs. 3 and 
4. 

The difference curve shows a maximum difference of 
less than 10 7 , which is far smaller than noise levels in real 
occultation data. However, it also shows two discontinui- 
ties. To investigate the origin of the discontinuities, we 
calculated the first and second differences of the flux from 
both the model and our method for the eight points sur- 
rounding each discontinuity. These differences showed that 
the second discontinuity appears in the flux from our 
method, but not in the flux from the model. The first dis- 
continuity did not seem to appear in either set of numbers, 
so we made a light curve from radii that were spaced 20 
times closer together around both discontinuities and found 
that the second one remained a discontinuity, but the first 



Scale heights from half light, (y-y h )/H 

Fig. 3 — Sample light curve. This light curve was calculated with the 
analytic benchmark model according to the prescription given in Appendix 
(A.l), for X* = 100 and r h = 1500. 


one clearly showed points spanning the gap of what had 
looked like a discontinuity in the original test. Both of these 
sharp changes are most likely due to numerical limitations 
in our method. 

3.1.2 Tightly Bound Atmosphere 

Since the first-order structure of an atmosphere is usually 
exponential rather than a power law, the second benchmark 
model is one in which the refraction angle is strictly an 
exponential function of radius (see Appendix A. 2). We ran 
a series of tests with this model. When we compared our 
method to the first-order approximation, V\(r) i we found 
that tests with X h = 4X10 7 showed the smallest differ- 
ences between our method and results calculated from the 
exponential refraction model, with larger differences at both 
larger and smaller values of X h (see Fig. 5). The differences 
for values larger than 4X 10 7 , shown in Fig. 6, were ran- 
domly scattered around zero, while differences for smaller 
values followed a smooth curve. This implies that the dif- 
ferences were due to round off at larger values and the 
first-order approximation at smaller values. 

When we re-ran the same set of tests with this model 
including the linear term, and again with this model includ- 
ing up through the cubic term, we got smaller errors but 
still found that there was a smooth difference curve at larger 



20 10 0 -10 -20 -30 -40 -50 -60 

Scale heights from half fight, (y-y^/H 

Fig. 4— Analytic test. Differences between our model results and the light 
curve in Fig. 3 (from the equations in Appendix (A.l) with X 0 = 100 and 
r h ~ 1500). The differences are defined to be flux from the model minus 
the flux from our method. Note the two apparent discontinuities in the 
difference curve — one around five scale heights below half Tight, and the 
other around 40 scale heights below half light. We believe that they are 
due to numerical limitations in the implementation of our method. 
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Fig. 5— Maximum size of differences in tightly bound atmosphere tests 
using different numbers of terms from the power series. Note that the 
differences in the first-order tests bottom out for \ h of about 4X10 , 
while the tests with more terms bottom out at slightly larger ratios and 
generally have smaller differences, especially at smaller values of \ h . The 
two sets of tests with more terms have similar differences, especially in the 
region where most real occulting bodies fall (X* of about 10 3x 10 ). 

ratios (see Fig. 7). In fact, this difference curve wus almost 
identical for these two sets of tests. 

The patterns we saw in the differences might be caused 
by problems with the power series. The refractivity along 
the path of a ray of light is Gaussian when x is small 
compared to r\ but farther out in the atmosphere where a: 
becomes comparable to r' t the path of the ray becomes 
similar to a radial path, and the refractivity approaches a 
simple exponential. In the case of a large planet, this change 
in the behavior of v occurs high in the atmosphere where v 
is very small and thus has little effect on the bending of the 
ray. In smaller planets, however, this change occurs lower 
in the atmosphere, and could cause the differences we saw 
in these tests. 



Fig. 6 — Differences in one term tightly bound atmosphere test with \ h 
= 4X 10 8 ( r h = 10 10 km, H = 25 km), defined as flux from model mi- 
nus flux from our method. The differences are well scattered with mean 
very close to zero. No pattern is visible. These differences are most likely 
caused by round off at the machine precision. 



Fig. 7 Differences in one term and two term tightly bound atmosphere 

tests with X* = 4x 10 3 (r h = 10 5 km, H = 25 km), as defined above. 
Deeper than about twelve scale heights below half light, ^ the two term 
differences dropped below zero, and were of order -10 11 . The differ- 
ences in both tests follow definite, smooth curves, with no evidence of 
scatter. However the curves have different shapes. Both of these shapes 
were consistent across \ h ratios, with the shape of the two term test also 
appearing consistently in tests where we included more than two terms. 

3.L3 Loosely Bound Atmosphere 

Next, we ran a series of tests using EY92*s isothermal 
model on sample planets similar to those in our solar sys- 
tem. Table 2 gives approximate half-light radii and scale 
heights for some of the planets and moons that have been 
observed by occultations. It shows the range of X values 
that naturally occur. We ran three groups of tests on model 
planets spanning the size range of the solar system, one 
similar to Jupiter, one similar to Venus, and one similar to 
Pluto. 

In each group of tests, we used model planets with four 
different combinations of values of half-light radius r h and 
energy ratio at half-light X* . We tested how the differences 
depended on several properties of the profiles: how many 
points per scale height we had in each profile, how far 
above the half-light radius the profiles extended, whether or 
not we extended the upper boundary with the isothermal 
model, and how many points were calculated for the inter- 
nal interpolation function for r(y). We did include the 


Table 2 

Half-Light Radii and Scale Heights for Planets and Satellites 
in the Solar System 3 


Planet or Satellite 

r h 

(km) 

H 

(km) 

x* 

Jupiter 

71900 

25 

2860 

Saturn 

61000 

62 

980 

Venus 

6200 

7 

890 

Neptune 

25300 

50 

500 

Uranus 

26100 

50 

480 

Mars 

3500 

8 

430 

Triton 

1450 

19 

77 

Titan 

3000 

47 

63 “ 

Pluto 

1200 

60 

20 r 


a Values were obtained from (Elliot and Olkin 1996) and references therein. 
The half-light radii are typical values for occulations observed from Earth, 
and scale heights are measured around the microbar pressure level in the 


atmosphere. 
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Table 3 

Maximum Differences between EY92\s Model and Our Method 3 



(a) Differing tops of refractivity profiles, upper boundary extended with 

an isothermal model 




Top of the model (x max ) in scale heights above half-light radius 


Model planet 

20 

15 

10 

5 

Jupiter-like 

8.5X10” 9 

2.5X10" 9 

4 X 10~ 8 

4 X 10 -6 

Venus-like 

1.5X10' 8 

8 X 10" 9 

I x 10 7 

9 X 10‘ 5 

Pluto-like b 

2X KT 5 

3 X 10“ 3 

3 X 10" 4 

2 X I0“ 3 


(b) Differing tops of refractivity profiles, upper boundary not extended with an isothermal model 




Top of the model (* max ) in scale heights above half flight radius 


Model planet 

20 

15 

10 

5 

Jupiter-like 

8.5X10' 9 

2.2x 10" 7 

3.5 x 10~ 5 

5.1 X 10“ 3 

Venus-like 

1 .3 X 1 0 * 8 

2.2X10 -7 

3.1X10' 5 

4.5XI0~ 1 

P!uto-like b 

I.9X 10“ ! 

6.2X 10’ 3 

4.5 X I0~ 4 

9.5X10' 3 


(c) Differing numbers of points per scale height in the refractivity profiles 




Points per scale height in refractivity profile 


Model planet 

10 

7 

5 

2 

Jupiter-like 

8.5X1CT 9 

8.5x 10~ 9 

8.5X I0 -9 

8.5X10" 9 

Venus-like 

I.3XI0" 8 

1.3X10" 8 

1.3X10" 8 

1 ,3X 10 8 

PIuto-like b 

2. OX 10 -3 

2.0x 10“ 5 

2.0X10’ 5 

I.9X10" 3 


(d) Differing numbers of points used in the internal interpolation of r{y) 




Points per scale height in refractivity profile 



Model planet 

125 

50 

10 


Jupiter-like 

8.5X 10 -9 

8.5X 10 9 

6.8X 10 -7 


Venus-like 

1.3X 10“ 8 

I.3X 10“ 8 

6.8 X 10" 7 


Pluto-like b 

l l j • . - r 

2.0 x 10 3 

2.0X 10 3 

2.1 x 10 -3 



3 Standard conditions for the tests were (i) to extend the profiles 20 scale heights above half light, (ii) to extend the upper boundary using the isothermal 
^atmospheric model, (iii) to have ten points per scale height in the profiles, and (iv) to calculate 125 points per scale height for r(y ) interpolation. 
Except for the last column, these numbers are averages of values that vary by a factor of five or more. 


focusing term in these tests and ran tests along nearly cen- 
tral chords. For simplicity, however, we did not test regions 
of the light curves near the central flashes. 

A summary of the results can be found in Table 3. All 
the results for the Jupiter- and Venus-sized model planets 



Fig. 8 — A typical difference plot for a Jupiter- or Venus-sized planet test 
vs. the loosely bound atmosphere model. The differences are concentrated 
in the region where the flux is dropping off the most rapidly. We extended 
the light curve out to 200 scale heights above half light and down to 450 
scale heights below' half light, but the differences were flat and very near 
zero in the regions not shown. This test was run with r h = 5000 km and 
H = 5 km (X A = 1 0 3 ) , and with a refractivity profile that had ten points 
per scale height and extended 20 scale heights above half light. We used 
125 points per scale height in calculating the interpolation function for 
r(y). 


(e.g., Fig. 8) were similar — none of the results were more 
than an order of magnitude apart. There was a weak trend 
for the differences in the Venus-sized cases to be larger 
than in the Jupiter-sized cases. In the Pluto-sized cases (e.g., 
Fig. 9), the differences were generally much larger and not 
as consistent from test to test as they were with the larger 
planets. 

In these tests, we consistently found that the most effi- 
cient form for the profile was to have only 2-5 points per 
scale height, but to have the profile extend about 15 scale 
heights above the half-light radius. Differences between our 
method and EY92’s model increased sharply when we low- 
ered the top of the profile (i.e., we decreased jc max ). This 
increase occurred more rapidly when we did not extend the 
upper boundary with the isothermal model than when we 
did. Minimum differences were achieved when we also 
calculated values for the interpolation of r(y) at 50 or more 
points per scale height — much finer spacing than in the 
profiles. 

The dependence of the size of the differences on the 
radius of the model planet may be explained by the power 
series approximation made in EY92’s model. In our Pluto- 
sized test planets, we used values of 10.48 and 20.97 for 
EY92’s model includes a power series in a variable 
proportional to 1 f\ h , and we found larger differences be- 
tween our method and their model in the cases where \ h 
was 10.48. The size of these differences could be explained 
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fio. 9 — A typical difference plot for a Pluto-sized planet test vs. the 
loosely bound atmosphere model. These differences are spread out over the 
entire region of the light curve, unlike the previous example. They also 
follow a much smoother curve. Unlike the larger planet tests, the differ- 
ences from each of these small-planet tests had their own characteristic 
error shape. These differences are most likely caused by approximations 
made in EY92’s model, not by errors in our method. This test ran with 
r A - 1250 km and H = 59.6 km (\ h = 20.8), and with a reactivity 
profile that had ten points per scale height and extended 20 scale heights 
above half light. We used 125 points per scale height in calculating the 
interpolation function for r(y). 

by the power series if the coefficient of the 5th term in their 
series is of order 10 . 

3.2 Multiple Light Curves 

In this section we discuss how our method works when 
we interpolate between profiles and when we calculate 
models within a least-squares fitting procedure. 

3.2, J Comparing Light Cur\>es 

We ran tests in which our method interpolated between 
profiles to calculate the light curve. We compared EY92’s 
no noise model for planets similar to Pluto and Triton 
(r h = 1500 km, X* = 100). For this type of planet, the 
errors for a single profile were around 10 7 . When we ran 
our method with three profiles spaced 10 % apart in \ h we 
got maximum differences between profiles of about 10 . 

When we shrank the spacing to 1% and 0.01%, the maxi- 
mum differences were about 10 7 in both cases the same 
as in the single profile case. When we ran similar tests with 
r h , we found that the maximum differences between pro- 
files were about 10 7 regardless of whether the profiles 
were spaced 10%, 1%, or 0.01% apart. This difference in 
behavior is most likely due to the linearity of r h — changing 
the half-light radius simply causes a linear shift in the light 
curve, while changing the energy ratio causes more compli- 
cated, nonlinear changes in the light curve. To verify this 
hypothesis, we compared light curves made with 10 % spac- 
ing in the square root, cube root, and natural logarithm of 
r h and got slightly larger differences than we did when we 
used r h itself as the parameter. 

3.2.2 Least-Squares Fits 

Next, we used EY92’s implementation of least-squares 
fitting in Mathematica™ to fit for \ h and r h in complete 


light curves with small central flashes. In fits for k h ancf r h 
with profiles spaced 1 % apart, the parameters converged to 
values that differed from the values used to generate the 
light curves by about seven parts in 10 7 errors on the 
same order as the differences between the model light 
curves and those generated by our method. The errors we 
found when we fit both parameters simultaneously were 
similar to the errors we found when we fixed one parameter 
at the model value and fit only for the other one. Fits for 
with profiles spaced 10 % apart yielded fitted values that 
differed from the model parameters by four parts in 10 - 
when the correct model parameter values fell between pro- 
files, and by seven parts in 10 7 when the model parameters 
were on a profile. Similar fits for r h yielded errors of about 
seven parts in 10 7 both on and off profile. 

To simulate fitting real data, we also ran a series of tests 
on EY92’s model with random Gaussian noise added to it. 
We characterized the noise level by specifying the ratio of 
the unocculted stellar signal to the rms noise level, inte- 
grated over one scale height of the atmosphere (measured in 
the shadow plane). We tested three different levels of noise 
that are typical of real data. The lowest signal-to-noise ratio 
per scale height used was 20 ; we also used values of 200 
and 2000 to simulate very high signal-to-noise ratios and to 
study how our method performed with data of different 
quality. In these tests, we spaced six profiles 3% apart in 
the one-parameter fits and 15 profiles (3 in r h by 5 in X/,) 
8 % apart in the two-parameter fits in order to have the 
profiles span the range of values the parameters converged 
to. A summary of our results is in Table 4. Note that 28 of 
42, or 67% of our results, fell within the formal error, 
which is very close to the 68 % expected. There appears to 
be a tendency for the r h fits to converge to a value smaller 
than the correct one — the fitted value was smaller than the 
model value in 15 of 21 tests. With this small a sample, 
however, the significance of this trend is difficult to deter- 
mine. Statistically, if there were no bias, there would be a 
4 % chance of getting results at least this skewed in this 
direction. 

As a comparison, we ran several fits of EY92 s model to 
noisy light curves generated the same way as we did in 
testing our method. These fits had very similar errors to 
those we got with our method. Upon fitting their model to 
light curves fit previously by our method, we found formal 
errors that were the same to within 1 % and fitted parameter 
values that were the same to within a couple percent of the 
formal errors. 

3.3 Computation Time 

To evaluate the speed of our method and verify that it is 
a faster alternative to calculating a complicated atmospheric 
model, we had Mathematica™ display the CPU time used 
in calculating the light curve each time we ran a test. For 
example for a Hewlett-Packard™ PA-RISC 7150 processor 
running at 125 MHz (SPEC CINT95 — 4.04, 
CF P95 = 4 . 55 ), the one-profile EY92 tests that we ran typi- 
cally took about 10 min of CPU time to calculate the light 
curve, ranging from 920 s down to 50 s. In these tests, the 
refractivity profiles extended 5-20 scale height above r h 
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Table 4 

Results of Fits to Noisy Data 


S/N = 20 


S/N = 200 


S/N = 2000 


Fit 

Fitted value 

Real 

value 

cr’s 

apart 3 

Fitted value 

Real 

value 

cr’s 

apart 8 

Fitted value 

Real 

value 

\ h alone 

103.4±6.6 

100.5 

+ 0.44 

I0I.92±0.62 

100.5 

+ 2.29 

100.626±0.065 

100.5 


104.1 ±6.8 

100.5 

+ 0.51 

99.08±0.61 

100.5 

-0.34 

I00.496±0.062 

100.5 


100.8±6.2 

101.2 

-0.07 

102.01 ±0.63 

101.2 

+ 1.26 

101.131 ±0.063 

101.2 


95.3±5.8 

99.7 

-0.77 

99.59±0.61 

99.7 

-0.18 

99.623 ±0.061 

99.7 

r h alone 

I501.0±2.6 

1507.5 

-2.51 

1507.11 ±0.26 

1507.5 

-1.51 

1507.51 2±0.026 

1507.5 


1504.9±2.6 

1507.5 

-0.98 

1507.00±0.25 

1507.5 

-2.01 

1507.488±0.026 

1507.5 


1519.9±2.5 

1521.0 

-0.42 

1520.82±0.26 

1521.0 

-0.68 

I520.976±0.023 

1521.0 


1497.2±2.6 

1498.0 

-0.32 

1498.09±0.25 

1498.0 

+ 0.38 

1498.003 ±0.025 

1498.0 

both, k h 

109.1 ±7.4 

100.5 

+ 1.17 

99.52±0.62 

100.5 

-1.59 

100.501 ±0.066 

100.5 


109.9±7.2 

101.2 

+ 1.20 

101 .20±0.65 

101.2 

0.00 

101. 159±0.063 

101.2 


99.7 ±6,7 

99.7 

0.00 

99.59±0.64 

99.7 

-0.17 

99.718±0.062 

99.7 

both, r h 

1508.1 ±2.5 

1507.5 

+ 0.20 

1507.74±0.24 

1507.5 

+ 1.03 

1507.491 ±0.026 

1507.5 


1519.1 ±2.4 

1521.0 

-0.79 

1520.70±0.26 

1521.0 

-1.18 

1520.992 ±0.025 

1521.0 


1494.4 ±2.7 

1498.0 

- 1.33 

1520.70±0.26 

1498.0 

-0.32 

1498.025±0.025 

1498.0 


“These numbers are (fitted value -real value)/(standard deviation of fitted value). 


cr’s 

apart 4 * * * 


+ 1.93 
-0.07 
- 1.09 
“1.27 
+ 0.48 
-0.46 
-0.96 
+ 0.10 
+ 0.01 
-0.64 
+ 0.29 
-0.36 
-0.33 
+ 0.98 


with 2-10 points per scale height, and we calculated 10- 
125 points per scale height in the internal interpolation 
function for r(y) and 200-700 points in the final light 
curve. The fastest tests were the ones where we did not 
calculate as many points for the internal interpolation of 
r(y). Extending the upper boundary with the isothermal 
model generally added about 25% to the calculation time. 
Using refractivity profiles with fewer points per scale height 
had little effect on the calculation time for the model, but 
lowering the top of the profile (decreasing x max ) increased 
the time considerably. The number of points at which we 
wished to calculate the light curve did not seem to be a 
major factor. Our tests often had different numbers of 
points in them, and some of the tests with more points in 
the light curve ran faster than tests with fewer. 

The time needed for least-squares fits depends on the 
number of parameters fit and the number of iterations 
needed for the values of the particular set of parameters to 
converge. In general, to determine the time needed for a fit, 
one should multiply the times given above by the number 
of parameters fit, and multiply that by the number of itera- 
tions needed. Our fits for two parameters often took as long 
as 3-4 h. If we implemented this method in C, we predict 
that such two parameter fits would take only about 2 min. 


4. CONCLUSIONS 

We have developed a method for numerically calculating 

occultation light curves from an arbitrary atmospheric 

model that is fast enough to enable least-squares fitting for 

model parameters. This method can be applied to a broad 
range of bodies, including Pluto and Triton. We have tested 
and verified that our method shows reasonable agreement 
with more exact calculations. Small, systematic errors in 
our method (most likely due to numerical procedures) are at 
a level much lower than the noise levels in recent occulta- 
tion observations. To date, the best signal-to-noise per scale 
height levels have been on the order of 10 3 . When we used 


sufficiently complete refractivity profiles, the largest differ- 
ences between our method and the comparison models were 
on the order of 1CT 4 . 

We found that it is better to interpolate the refractivity 
derivatives than to interpolate light curves, and that the 
most efficient form for refractivity profiles that still gave 
good results was for each profile to extend well above the 
half-light radius but to include only a few points per scale 
height. The quality of results obtained with such profiles 
implies that, if necessary, it is acceptable to record occul- 
tation data with as few as two light-curve points per scale 
height. 

Planned improvements to our procedure include speeding 
up the calculations by implementing the code in a compiled 
language, such as C. This will enable the least-squares fit- 
ting of atmospheric models with more parameters. Also, for 
large, oblate planets — which do not satisfy the assumption 
of spherical symmetry — one could recast the fundamental 
equations in terms of local radii of curvature, rather than the 
single radius used for the spherical-symmetry assumption. 

We thank Catherine Olkin for interesting discussions, 
helping us debug our code, and providing incentive for us 
to get it working quickly when she needed it for her thesis. 
We also are grateful to Bill Hubbard for his comments as 
referee and appreciate suggestions from Mary Agner. This 
work was supported, in part, by NASA Grant NAGW-1494, 
NSF Grant AST-9322115, and NSF’s Research Experiences 
for Undergraduates (REU) Program. 


APPENDIX: EQUATIONS USED FOR MODEL- 
ATMOSPHERE TESTS 

In this Appendix we develop the equations in the form 
needed for the tests of our method performed on an analytic 
atmosphere and a tightly bound atmosphere. Equations for 
the loosely bound atmosphere appear in EY92. 
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A.l Analytic Atmosphere 

Eshleman and Gurrola (1993) have presented a set of 
descriptive atmospheric models that satisfy hydrostatic equi- 
librium for different power-law radial dependencies of grav- 
ity and temperature. For a spherically symmetric gravity 
field (i.e., g(r) = GM/r 2 ), an inverse radial temperature 
dependence (i.e., T(r) = T 0 r 0 /r , where r 0 is a reference 
radius) yields a model expressible with ordinary functions. 

If H 0 is the pressure scale height at the reference radius, 
we define X 0 = r 0 /H 0 . The refractivity and its first two 
derivatives for this analytic model, needed to calculate the 
light curve with our method, are given by 

I r\^ o-’) 

Hr) = 

dv{r) _ Vq(Xq-I) / r _ \ K ° 
dr r 0 \r 0 ) 


(Al.l) 
(A 1.2) 


d 2 HH 

7? 


vo^oi^o 1 ) l r _\ ! 0 

Z2 rj 


(A 1.3) 


We also need the analytic expressions for the refraction 
angle and its derivative, required for calculating an analytic 
light curve (Eshleman and Gurrola 1993): 


6{r) = - 


2 v 0 TT m r\ T / r\-^o-D 


r ~2 2 


Its derivative is 


d6(r) 


r o 


^•o 


dr 


2n 0 7T 1/2 (X 0 — 1 ) r ( — 


Xn 1 

r » r ly -2 


(A 1.4) 


r 

r ol 


(A 1.5) 


The value of the reference radius, r 0 , is arbitrary, so we 
choose to set it to the half-light radius, r h . To find r h , we 
set 7 cyi to 0.5 in Eq. (4) and solve for dOldr. Substituting 
this expression into Eq. (A 1.5) and solving for r 0 gives us 
the following formula for r h : 


r o = r h ~ 


2 v 0 tt v 2 d(\ 0 - or 


r) 



(A 1.6) 


then, by definition, X at half-light, X/, = Xq- To calculate 
the analytic light curve, we substitute Eq. (A1.6) for r 0 in 
Eq. (A 1.5) and then put this equation into Eq. (4). Equa- 
tions (5) and (A 1.4) can be used to relate r and y. 


A.2 Tightly Bound Atmosphere 

In the exact equation for a cylindrical atmosphere (one 
with no focusing in the plane perpendicular to the ray path), 
the light-curve flux, 7cy\ , is related to the derivative of the 
refraction angle, 9(r) for an atmosphere at a distance D by 
Eq. (4): 


7cy\ 


] ■ n 


do 

77 


(A2.1) 


Also, y is related to the refraction angle in the small-angle 
approximation by Eq. (5): 

y(r) = r+D9(r). (A2.2) 

If we pick a reference ray that intercepts the atmosphere at 
a radius r 0 and is then refracted by the atmosphere by an 
angle 9 0 , so that it intercepts the observer at y 0 , we have 

y 0 =r 0 +D9 0 . (A2.3) 

We postulate that the atmosphere produces an exponential 
refraction angle with a scale height H, so that 

9(r) = 9 0 e- (r ~ r ° )IH . (A2.4) 

Next, we subtract (A3) from (A2) and divide by the scale 
height H : 

(y-yo) = l { o-e 0 ). (A2.5) 

H ti ti 


We choose to put the reference level at half light 
(<£ cy | = 1/2). Now we can solve for 0 O : 


d* cy] 'y 


de 

l+D — I 
dr 


1 

~d¥o' 

H 

(A2.6) 


or 


H 

00 = ~ D' 


(A2.7) 


Substituting 9 0 = - HID in to the flux equation (Eq. (4)), 
we get 

i 1 _ 1 

& cyi = 1 Jo\ ~ |l +e“ (r-r ° )/H | 1 +e _(r_r ° ,/H 


d9 

1 + 0 37 1 


(A2.8) 

Solving for the exponential in terms of the flux, we find 


e -{r-r 0 )!H - _i J 

4> cyl 


The exponent is 


(r~r 0 ) 

H 


= -In 


( 7 cyl ) ’ 


(A2.9) 


(A2.10) 


which is the first term on the right-hand side of Eq. (A2.5). 
Next we look at the second term on the right-hand side of 
Eq. (A2.5): 

D (f) _0) ~ — — = 111 r e -(r-r 0 )/H_,] 

H (® 6 °) " H H H L 


7 cyl 


(A2.ll) 
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where the last step in Eq. (A2.ll) used Eqs. (A2.4) and 
(A2.9). Substituting this result and Eq, (A2.10) into Eq. 
(A2.5) and denoting this is the observer’s y coordinate at 
half light by y h , we find 

y(0-yh = I 1 ,\-f — 

H Ucyl(y) J Ucyl (y) 

(A2.12) 

This light-curve equation, valid for exact exponential refrac- 
tion, is identical to the approximation derived by Baum and 
Code (1953) for the occultation light curve for a large, 
spherical planet with an isothermal atmosphere. 



Next, we need a formula for the refractivity, v(r), that 
yields an exponential 0(r) given by Eq. (A2.4). To find the 
desired function we first tried a simple exponential and 
calculated the refraction angle 6 for this v. We found that 0 
would be exponential to first order if we multiplied the 
exponential by 1 fyfr: 

v x {r) = v h yj^- e - (r - r * )m . (A2.13) 

To make 0 be exponential to higher order, we multiplied 
v ! by a power series in Hfr , choosing the coefficients to 
cancel out the nonexponential terms in 0(r): 


Hr) = v h 



J - 


1 

8 



H\ 2 75 (H\ 3 3675 ///\ 4 59535 I H\ 5 

7 / ~~ l024 \ 7 / + 32768 \ 7 / ” 262144 \ 7 / + "' ' 

(A2.14) 


Note that the coefficients begin to slowly increase in 
size, indicating that this series is asymptotic and will most 
likely diverge if we include enough terms. 

REFERENCES 

Baum, W. A., and Code, A. D. 1953, AJ, 58, 108 
Born, M., and Wolf, E. 1964, Principles of Optics (New York, 
MacMillan) 

Chamberlain, D. M. 1996, A Numerical Method for Calculating 
Occultation Light Curves from an Arbitrary Atmospheric 
Model, S.B. thesis, Massachusetts Institute of Technology 
Elliot, J. L., et al. 1993, AJ, 106, 2544 

Elliot, J. L., French, R. G., Dunham, E., Gierasch, P. J., Veverka, 
J., Church, C., and Sagan, C. 1977, ApJ, 217, 661 
Elliot, J. L., and Olkin, C. B. 1996, in Annual Review of Earth 
and Planetary Sciences, ed. G. W. Wetherill (Palo Alto, Annual 


Reviews), p. 89 

Elliot, J. L., and Young, L. A. 1992, AJ, 103, 991 
Eshleman, V. R., and Gurrola, E. M. 1993, Icarus, 105, 298 
French, R. G., Elliot, J. L., and Gierasch, P. J. 1978, Icarus, 33, 
186 

Kovalevsky, J., and Link, F. 1969, A&A, 2, 398 
Krasnopolsky, V. A,, Sandel, B. R., Herbert, F., and Vervack, R. 
J. 1993, JGR, 98, 3065 

Olkin, C. B. 1996, Ph.D. thesis, Massachusetts Institute of Tech- 
nology 

Roques, F., et al. 1994, A&A, 288, 985 

Strobel, D. F., and Summers, M. E. 1995, in Neptune and Triton, 
ed. D. P. Cruikshank (Tucson, University of Arizona Press), p. 
1107 

Wasserman, L. H,, and Veverka, J. 1973, Icarus, 20, 322 
Wolfram, S. 1991, Mathematica (Redwood City, CA, Addison- 
Wesley) 



